Transverse metacenter height estimation device and transverse metacenter height estimation method

ABSTRACT

Provided is a transverse metacenter height estimation device that includes: a history storage means for storing time sequence data of the roll angle of a hull; and a transverse metacenter height estimation means that estimates the transverse metacenter height of the hull on the basis of the time sequence data of the roll angle of the hull, such data stored by the history storing means. The transverse metacenter height estimation device is characterized in that the transverse metacenter height estimation means first calculates a roll natural frequency on the basis of the time sequence data of the roll angle of the hull, sets the calculated roll natural frequency as an observation model, carries out state estimation on the basis of a normal state spatial model in which the transverse metacenter height and the gyration radius of the hull are used as state variables, and thereby estimates the transverse metacenter height.

INCORPORATION BY REFERENCE

The present application claims priority from Japanese Patent ApplicationNo. 2014-104786 filed on May 20, 2014, the contents of which are herebyincorporated by reference.

TECHNICAL FIELD

The present invention relates to a transverse metacenter heightestimation device and a transverse metacenter height estimation method.

BACKGROUND ART

Conventionally, with respect to a vessel running on randomly movingwaves, it is important to appropriately grasp hull motion data, hullstate data, and sea conditions from a safety point of view. The hullmotion data are the data related to the motion of the hull, such as adisplacement and an acceleration of the hull. Further, the hull statedata are the data related to the state of the hull, such as a draft, adrainage volume, and a transverse metacenter height (hereinafter, alsoreferred to as “GM”) of the hull. Further, the sea conditions are theinformation related to the sea phenomena, such as a wave height, a wavecycle, and a wave direction of the wave on the area where the vessel isrunning.

Conventionally, it is common that the hull motion data of these piecesof information are kinetically obtained on the basis of the previouslyset hull state data and the sea conditions provided from an informationproviding institution such as the Meteorological Agency.

However, in this method, it is impossible to grasp the hull motion dataand the hull state data appropriately. This is because the seaconditions provided from an information providing institution such asthe Meteorological Agency have a large amount of information on a widearea of water instead of a local area of water in which the vessel iscurrently running, and have a low accuracy.

To address this issue, in recent years, a method is reported in whichhull motion data are measured with various devices mounted on a vessel,the measured hull motion data, which are unsteady time-series data, aresubjected to statistic processing in real time, so that the hull statedata and the sea conditions are statistically estimated (for example,see Non-Patent Document 1).

Non-Patent Document 1 discloses a technique in which the sea conditionsare estimated by analyzing the hull motion data, which are unsteadytime-series data, by using the Time Varying Coefficient Vector AR(TVVAR).

CITATION LIST Non-Patent Literature

Non-Patent Literature 1: Tsugikiyo Hirayama, Toshio Iseki, ShigesukeIshida, “Real-Time Detection Method of Encountering Ocean Waves andRecent Results”, The Japan Society of Naval Architects and OceanEngineers, December 2003, pp. 74-96.

SUMMARY OF THE INVENTION Technical Problems

As described in the above conventional art, the hull state data areestimated on the basis of the hull motion data as the unsteadytime-series data. As an example, “GM”, which is one of the hull statedata, is estimated on the basis of “time-series data of a rollingangle”, which is one of the hull motion data. Specifically, the rollingnatural frequency is first estimated on the basis of the time-seriesdata of the rolling angle, and the GM is then estimated on the basis ofthe estimated rolling natural frequency.

However, by the conventional estimation method, it is impossible toaccurately estimate the GM. This is because when the GM is estimated onthe basis of the rolling natural frequency, an approximation formula asshown in Equation (1) is used.

$\begin{matrix}\left\lbrack {{Mathematical}\mspace{14mu} {Expression}\mspace{14mu} 1} \right\rbrack & \; \\{T = \frac{2\; {CB}}{\sqrt{GM}}} & (1) \\{C = {0.373 + {0.023\left( \frac{B}{d} \right)} - {0.043\left( \frac{L}{100} \right)}}} & \;\end{matrix}$

In Equation (1), T is a rolling natural period, C is an experimentalconstant, B is a width of a ship, d is a draft of the ship, and L is alength of the ship. In Equation (1), C is used as a constant, where Cindicates a different value depending on loading conditions and thelike. For this reason, the accuracy of the estimated value of GM isaccordingly low.

The present invention has been made in view of the above issue, and anobject of the invention is to provide a transverse metacenter heightestimation device and a transverse metacenter height estimation methodwhich enable highly accurate estimation of the transverse metacenterheight.

Solution to Problems

In order to achieve the above object, an transverse metacenter heightestimation device according to the present invention includes: a timehistory memory unit that stores time-series data of a rolling angle of ahull; and a transverse metacenter height estimation unit that estimatesa transverse metacenter height of the hull based on the time-series dataof the rolling angle of the hull stored in the time history memory unit.The transverse metacenter height estimation unit first calculates arolling natural frequency based on the time-series data of the rollingangle of the hull, and then estimates, while using the calculatedrolling natural frequency as an observation model, the transversemetacenter height by performing state estimation based on a generalstate-space model in which the transverse metacenter height and a radiusof gyration of the hull are state variables.

Advantageous Effects of Invention

The present invention enables estimation of a transverse metacenterheight with high accuracy.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a diagram showing an example of a hardware configuration of asea phenomena estimation system including a transverse metacenter heightestimation device according to an embodiment.

FIG. 2 is a diagram showing a functional configuration example of thesea phenomena estimation system including a transverse metacenter heightestimation device according to the embodiment.

FIG. 3 is a flowchart showing an example of a control logic of the seaphenomena estimation system including a transverse metacenter heightestimation device according to the embodiment.

DESCRIPTION OF EMBODIMENT

Hereinafter, an embodiment according to the present invention will bedescribed.

Hardware Configuration of System

FIG. 1 is a diagram showing an example of a hardware configuration of asea phenomena estimation system including a transverse metacenter heightestimation device according to the embodiment.

A sea phenomena estimation system 1 shown in FIG. 1 is equipped with asatellite compass 2, an information processing device (transversemetacenter height estimation device) 3, and a display 4. The seaphenomena estimation system 1 is mounted in a hull of a vessel.

The satellite compass (GPS compass) 2 is a device having a function as adirection sensor which calculates the direction of the hull on the basisof the relationship between relative positions of two GPS antennasattached in a bow direction of the hull. The satellite compass 2 hasalso a function as a motion in wave sensor which can measure transversemotion in wave (rolling), longitudinal motion in wave (pitching), andvertical motion in wave (heaving) of the hull. Note that a gyro sensormay be used instead of the satellite compass 2.

The information processing device 3 is a computer device which includesa memory device 31, a processing unit 32, an interface device 33, aninput device 34, an auxiliary storage device 35, and a drive device 36which are each connected to one another through a bus 38. Theinformation processing device 3 estimates the transverse metacenterheight on the basis of information measured by the satellite compass 2.Further, the information processing device 3 estimates the seaconditions on the basis of the estimated value of the transversemetacenter height and the like. The information processing device 3corresponds to a “transverse metacenter height estimation device” of theclaims. Note that the information processing device 3 and the display 4to be described later may be integrated with the satellite compass 2.

The memory device 31 is a storage device such as a random access memory(RAM) which reads out, at a time of start-up of the informationprocessing device 3, a program (a program which realizes a hull statedata calculation section 23 and a sea phenomena estimation section 24 ofFIG. 2) and the like stored in the auxiliary storage device 35. Thememory device 31 stores also a file, data, and the like necessary toexecute a program.

The processing unit 32 is a processing unit such as a central processingunit (CPU) which executes a program stored in the memory device 31. Theinterface device 33 is an interface device to connect to an externaldevice such as the satellite compass 2 and the display 4. The inputdevice 34 is an input device (for example, a keyboard or a mouse) toprovide a user interface.

The auxiliary storage device 35 is a storage device such as a hard diskdrive (HDD) which stores a program, a file, data, and the like. Theauxiliary storage device 35 stores a program and the like which realizea function of the hull state data calculation section 23 and the seaphenomena estimation section 24 of FIG. 2.

The drive device 36 is a device which reads out a program (for example,a program which realizes a function of the hull state data calculationsection 23 and the sea phenomena estimation section 24 of FIG. 2) storedin a storage medium 37. A program read out by the drive device 36 isinstalled in the auxiliary storage device 35. The storage medium 37 is astorage medium such as a universal serial bus (USB) memory and an SDmemory card, and the storage medium store the above program or the like.

The display 4 is an output device which outputs, on a screen, outputdata generated by the information processing device 3, for example, thesea conditions.

Functional Configuration of System

FIG. 2 is a diagram showing a functional configuration example of thesea phenomena estimation system including a transverse metacenter heightestimation device according to the present embodiment. Note that, in thefollowing description, components similar to the components of FIG. 1are assigned the same reference numerals, and redundant description isappropriately omitted.

The sea phenomena estimation system 1 shown in FIG. 2 has a measurementsection 21, a time history memory 22, a hull state data calculationsection (transverse metacenter height estimation section) 23, the seaphenomena estimation section 24, and an output section 25. The seaphenomena estimation system 1 is mounted on a vessel. Note that theinformation processing device 3 realizes respective functions of thetime history memory 22, the hull state data calculation section 23, andthe sea phenomena estimation section 24.

The measurement section 21 is a measurement unit which measures hullmotion data of a vessel on which the sea phenomena estimation system 1is mounted. Here, the hull motion data means the data related to themotion of the hull such as a rolling angle, a pitching angle, adisplacement of heaving, and the like of the hull. Note that it is alsopossible to use angular velocities of rolling and pitching and anacceleration of heaving of the hull. The measurement section 21 isrealized by the satellite compass 2 of FIG. 1 or a gyro sensor.

The time history memory 22 is a time history storage unit which stores atime history of the hull motion data measured by the measurement section21. The time history memory 22 stores time-series data of the hullmotion data in a predetermined period from past to now. The time historymemory 22 is realized by the memory device 31 and the like of FIG. 1.Note that an input source of the history of the hull motion data is notlimited to the measurement section 21. For example, it is also possiblethat another information processing device storing a history of the hullmotion data is used as the input source.

The hull state data calculation section 23 is a hull state datacalculation unit which calculates the hull state data on the basis ofthe history of the hull motion data stored in the time history memory 22in a predetermined period. Here, the hull state data means the datarelated to the state of the hull such as a draft, a drainage volume, anda GM. The hull state data calculation section 23 is realized by theprocessing unit 32 and the like of FIG. 1. The hull state datacalculation section 23 corresponds to a “transverse metacenter heightestimation unit” of the claims.

The sea phenomena estimation section 24 is a sea phenomena estimationunit which estimates local sea conditions in the water in which thevessel equipped with the sea phenomena estimation system 1 is running,on the basis of the history of the hull motion data stored in the timehistory memory 22 and the hull state data calculated by the hull statedata calculation section 23. Here, the sea conditions means theinformation related to the sea phenomena such as a wave height, a wavecycle, a wave direction, and the like of the wave in the area in whichthe vessel is running.

The sea phenomena estimation section 24 is realized by the processingunit 32 and the like of FIG. 1.

The output section 25 is an output unit which outputs the hull statedata calculated by the hull state data calculation section 23 and thesea conditions estimated by the sea phenomena estimation section 24. Theoutput section 25 is realized by the display 4 and the like of FIG. 1.

With the configuration described above, in the sea phenomena estimationsystem 1 according to the present embodiment, the hull state datacalculation section (transverse metacenter height estimation section) 23owned by the information processing device 3 estimates the transversemetacenter height on the basis of the hull motion data measured by themeasurement section 21, and then the sea phenomena estimation section 24estimates the sea conditions. An output section 25 outputs the estimatedtransverse metacenter height and the estimated sea conditions.

Control Logic of System

FIG. 3 is a flowchart showing a control logic of the sea phenomenaestimation system including a transverse metacenter height estimationdevice according to the present embodiment.

The sea phenomena estimation system 1 repeatedly estimates the seaconditions by repeatedly performing a control logic of a series of stepsS1 to S8 shown in FIG. 3. In particular, the sea phenomena estimationsystem 1 estimates the transverse metacenter height by repeatedlyperforming the process of steps S1 to S3. Note that a description willbe given below, appropriately with reference to FIG. 2.

First, in step S1, the measurement section 21 measures the hull motiondata (step S1). Specifically, the data of the rolling angle and thepitching angle and the displacement of heaving of the hull. Step S1 maybe successively performed in the process of repeating the series ofsteps S1 to S8, or may be repeatedly performed, by a batch process orthe like, as a process independent from steps S2 to S8.

Note that by repeatedly performing the process of step S1 shown in FIG.3, the time-series data of the rolling angle, the pitching angle, andthe displacement of heaving in a predetermined period from past to now,in other words, the data of the angles and the like at respective timesare stored in the time history memory 22. Then, the process goes torespective processes of steps S2, S4, and S6.

Next, in step S2, the hull state data calculation section 23 calculates(estimates) a rolling natural frequency on the basis of the time-seriesdata of the rolling angle (step S2). Although the process of step S2 isa known technique, an example of the process will be described below.

Specifically, a second order linear probability dynamic model (seeEquation (2) below) about time-series data x(t) of rolling isconsidered. Note that, in Equation (2), a₁(=2α) is a damping coefficientand a₂(=ω²) is a square of a natural angular frequency ω. The term u(t)represents an external force term dealt as a stochastic process and hasa finite dispersion. However, the term u(t) does not need whiteness.

[Mathematical Expression 2]

x″(t)+a ₁ x′(t)+a ₂ x(t)=u(t)   (2)

Further, the external force term u(t) in Equation (2) is expressed by anm order continuous auto regression model shown in Equation (3) below.Note that, in Equation (3), b_(i)(i=1, . . . , m) is a coefficient ofthe model, and v(t) is a normal white noise, where the average is 0, andthe dispersion is σ².

$\begin{matrix}\left\lbrack {{Mathematical}\mspace{14mu} {Expression}\mspace{14mu} 3} \right\rbrack & \; \\{{{u^{(m)}(t)} + {\sum\limits_{i = 1}^{m}{b_{i}{u^{({m - i})}(t)}}}} = {v(t)}} & (3)\end{matrix}$

By substituting Equation (2) into Equation (3), a whitened (m+2) ordercontinuous auto regression model shown in Equation (4) is obtained. Notethat c_(i) (where i=1, . . . , m+2) in Equation (4) is a coefficient ofthe model.

$\begin{matrix}\left\lbrack {{Mathematical}\mspace{14mu} {Expression}\mspace{14mu} 4} \right\rbrack & \; \\{{{x^{({m + 2})}(t)} + {\sum\limits_{i = 1}^{m + 2}\; {c_{i}{u^{({m + 2 - i})}(t)}}}} = {v(t)}} & (4)\end{matrix}$

Equation (4) is expressed as Equation (5) below in a vector form.

$\begin{matrix}\left\lbrack {{Mathematical}\mspace{14mu} {Expression}\mspace{14mu} 5} \right\rbrack & \; \\{{{x^{\prime}(t)} = {{{Ax}(t)} + {{Bv}(t)}}}{where}\left\{ \begin{matrix}{{x(t)} = \left( \; \begin{matrix}{x(t)} & {x^{\prime}(t)} & {x^{''}(t)} & \ldots & x^{({m + 1})}\end{matrix} \right)^{T}} \\{A = \begin{pmatrix}0 & 1 & 0 & \ldots & 0 \\0 & 0 & 1 & 0 & \vdots \\\vdots & \vdots & \ddots & \ddots & 0 \\0 & 0 & \ldots & 0 & 1 \\{- C_{m + 2}} & {- C_{m + 1}} & \ldots & {- C_{2}} & {- C_{1}}\end{pmatrix}} \\{{B = I_{{({m + 2})} \times {({m + 2})}}},} \\{{v(t)} = \begin{pmatrix}0 & 0 & \ldots & 0 & {v(t)}\end{pmatrix}^{T}}\end{matrix} \right.} & (5)\end{matrix}$

Equation (5) is dealt as a system model of a state space model afterbeing discretized. Further, in order to simultaneously estimate theunknown coefficients c_(i) (where i=1, . . . , m+2), a state vector ofthe state space model is considered with the above-described unknowncoefficients a_(i) and b_(i) included therein, and the model is thenextended into a auto-organizing state space model; then, the stateestimation and the unknown coefficients are simultaneously estimated byusing an Ensemble Kalman Filter. The state estimation by using theEnsemble Kalman Filter is a known technique and is thus not describedhere.

By the procedure described above, in step S2, the rolling naturalfrequency is calculated (estimated) on the basis of the time-series dataof the rolling angle stored in the time history memory 22. Note that therolling natural frequency may be calculated by a method different fromthe above method. For example, the rolling natural frequency may becalculated by using a discretized auto regression model.

After that, in step S3, the hull state data calculation section 23calculates (estimates) the GM on the basis of the rolling naturalfrequency calculated in step S2 (step S3).

In step S3, in a non-linear observation model in which the rollingnatural frequency (or a rolling natural period, which is the inverse)calculated in step S2 is used as observed data and in which the GM andthe radius of gyration are the state variables, it is assumed that thestate variables fluctuate slightly with time, and this is considered asthe system model; thus, a general state-space model analysis isperformed to simultaneously estimate the GM and the radius of gyration.

That is, the relationship represented by Equation (6) below holdsbetween the rolling natural frequency f (=ω/2π) calculated in step S2and the GM. Note that, in Equation (6), T is the rolling natural period(unit: s), f is the rolling natural frequency (unit: Hz), k is theradius of gyration (unit: m), and g is the gravitational acceleration(unit: m/s²). In Equation (6), GM (unit: m) and k are unknown.

$\begin{matrix}\left\lbrack {{Mathematical}\mspace{14mu} {Expression}\mspace{14mu} 6} \right\rbrack & \; \\{T = {\frac{1}{f} = \frac{2\pi \; k}{\sqrt{gGM}}}} & (6)\end{matrix}$

Then, in step S3, the state space model is considered in which twounknowns of the GM and the radius of gyration k are the state variables.That is, a non-linear observation model is considered in which therolling natural period Tn is observed from the estimated amount of thestates of GMn and kn at time n. It is assumed that the state variablesfluctuate slightly with time, and this is considered as the systemmodel. Thus, the general state-space model represented by Equation (7)below is configured. Note that, in Equation (7), vn is a system noise,and wn is an observation noise. For the sake of simplicity, the noisesare considered as normal white noises.

$\begin{matrix}\left\lbrack {{Mathematical}\mspace{14mu} {Expression}\mspace{14mu} 7} \right\rbrack & \; \\\left\{ {\begin{matrix}{x_{n} = {x_{n - 1} + v_{n}}} \\{T_{n} = {h\left( {x_{n},w_{n}} \right)}}\end{matrix}{where}\left\{ \begin{matrix}{x_{n} = \left( {{GM}_{n}\mspace{14mu} k_{n}} \right)^{T}} \\{v_{n} = \left( {v_{1,n}\mspace{14mu} v_{2,n}} \right)^{T}}\end{matrix} \right.} \right. & (7)\end{matrix}$

In step S3, state estimation is performed by using a Monte Carlo filter,which is a type of particle filters, on the basis of the generalstate-space model represented by Equation (7). Note that the stateestimation by using a Monte Carlo filter is a known method and is notdescribed here.

By the process of steps S2 and S3 described above, the hull state datacalculation section 23 calculates (estimates) the GM of the hull in realtime by analyzing the time-series data of the rolling angle in apredetermined period stored in the time history memory 22. By thismethod, it is possible to grasp the change in the position of the centerof gravity in a loading state and grasp a degree of motion of thevessel.

Further, because no approximation is used to estimate GM, it is possibleto obtain a generally stable estimation result of GM. For example, inthe case that GM of a vessel at the time of design was 0.52, GM wascalculated to be about 0.68 by the conventional method represented byEquation (1), but GM was estimated to be in the range of 0.48 to 0.54 bythe method according to the present embodiment. Therefore, the methodaccording to the present embodiment enables GM to be estimated with highaccuracy. Note that in the above Equation (7), the equation may be madeto include a valuable older than the previous timing (for example,X_(n-2)).

Returning to FIG. 3, if the process goes from step S1 to step S4, thehull state data calculation section 23 calculates the draft and thedrainage volume of the hull by analyzing the time-series data of thedisplacement of heaving stored in the time history memory 22 in step S1(step S4). For example, the draft is calculated on the basis of thetime-series data of the displacement of heaving, an installation heightof the GPS antennas (corresponding to the measurement section 21 of FIG.1), an inclination angle of the vessel in the longitudinal direction,and the like. Then, the drainage volume is further calculated on thebasis of the calculated draft.

If the process goes from steps S3 and S4 to step S5, the sea phenomenaestimation section 24 calculates a current hull response of the hull onthe basis of the hull state data calculated in step S2 (step S5).

As the process of step S5, there are two possible methods describedbelow, and any one of them is used. In the first method, a hull responsefunction is previously calculated to make a data base by using asparameters the hull state data (the draft and GM), a vessel speed, andsea phenomena (a wave height, a wave cycle, and a wave direction) whichis to be an input; and the hull response function corresponding to thecurrent states is then obtained by an interpolation calculation. In thesecond method, the hull response function is calculated on the basis ofa calculation formula, with the current state of the hull and the seaphenomena being used as inputs. In any of the methods, the sea phenomenaare an unknown to be obtained below and are also a term necessary tocalculate the hull response function. The most appropriate responsefunction is selected in real time as a non-linear problem by using aniteration method. Note that the response function is here a functionindicating how the hull responds (moves) when the hull receives waveswith a regular wavelength from any given direction, and the parametersof the function are a wave direction, a wave length, and the like.

In the process of step S5, the most appropriate hull response functionis selected in real time, depending on the current hull motions (thepitching, the rolling, and the heaving) and the vessel speed measured bythe measurement section 21 (the satellite compass 2). By this process,even on the actual sea in which the vessel speed changes from moment tomoment, the most appropriate response function of the hull can beobtained. As a result, it is also possible to improve the accuracy ofthe estimation of the sea conditions to be described later.

Further, if the process goes from step S1 to step S6, the hull statedata calculation section 23 calculates cross spectra of the respectivehull motions (the rolling, the pitching, and the heaving) on the basisof the displacement of heaving, the pitching angle, and the rollingangle stored in the time history memory 22 in step S1 (step S6). Theprocess of step S6 is a known technique and is thus not described here.

In step S6, on the basis of the time-series data of the displacement ofheaving (unit: m), the pitching angle (unit: rad), the rolling angle(unit: rad), the hull state data calculation section 23 calculates foreach frequency, cross spectra of the respective hull motions (therolling, the pitching, and the heaving) including a rolling autospectrum (unit: rad²/s), a pitching auto spectrum (unit: rad²/s), anheaving auto spectrum (unit: m²/s), a pitching-heaving cross spectrum(unit: rad·m/s), a rolling-heaving cross spectrum (unit: rad·m/s), and apitching-rolling cross spectrum (unit: rad²/s). The cross spectraobtained for respective frequencies are stored as the time-series datain the time history memory 22.

If the process goes from steps S5 and S6 to step S7, the sea phenomenaestimation section 24 probability statistically calculates a directionalwave spectrum on the basis of the current hull response calculated instep S5 and the cross spectra of the respective hull motions (therolling, the pitching, and the heaving) calculated in step S6 (step S7).Hereinafter, such process will be described together with a theory, usedin step S7, according to the present embodiment.

If it is assumed that an ocean wave is represented by superposingcomponent waves coming from every direction and having all frequencies,a sea surface variation amount η(t) at time t at an fixed point (vesselposition) is expressed by Equation (8) below using a directional wavespectrum E(f, x) (unit: m²/(rad/s)). Note that, in Equation (8), thepart under the root symbol and ε(f, x) are respectively the amplitudeand the phase of a component wave having a frequency f coming from adirection x.

[Mathematical Expression 8]

η(t)=∫_(-π) ^(π)∫₀ ^(on) cos{2πft+ε(f,x)}√{square root over (2E(f,x)dfdx)}  (8)

On the other hand, if it is assumed that the hull motion in waveresponds linearly to an input wave, the directional wave spectrumE(f_(e), x) at an encounter frequency f_(e) of one wave, and a crossspectrum φln(f_(e)) of the hull motion in wave is express generally byEquation (9) below. Note that, in Equation (9), l and n are the modes ofthe hull motion in wave, and H_(l)(f_(e), x) and H_(n)*(f_(e), x) arerespectively the response functions in the modes l and n. Further, x isan encounter angle with respect to a wave, and the symbol “*” representsa complex conjugate.

[Mathematical Expression 9]

Φ_(in)(f _(e))=∫_(-π) ^(π) H _(l)(f _(e) , x)H*_(n)(f _(e) , x)E(f _(e), x)dx   (9)

Because Equation (9) is expressed based on the encounter frequency, thisequation is converted into an equation based on an absolute frequency(see Equation (10) below).

$\begin{matrix}{\mspace{79mu} \left\lbrack {{Mathematical}\mspace{14mu} {Expression}\mspace{14mu} 10} \right\rbrack} & \; \\{{\varphi_{\ln}\left( f_{e} \right)} = {{\int_{\frac{\pi}{2}}^{\pi}{{H_{l}\left( {f_{01},x} \right)}{H_{n}^{*}\left( {f_{01},x} \right)}{E\left( {f_{01},x} \right)}{\frac{f_{01}}{f_{e}}}\ {x}}} + {\int_{- \frac{\pi}{2}}^{\frac{\pi}{2}}{{H_{l}\left( {f_{01},x} \right)}{H_{n}^{*}\left( {f_{01},x} \right)}{E\left( {f_{01},x} \right)}{\frac{f_{01}}{f_{e}}}\ {{x\left( {f_{e} < \frac{1}{4\; A}} \right)}}}} + {\int_{- \frac{\pi}{2}}^{\frac{\pi}{2}}{{H_{l}\left( {f_{02},x} \right)}{H_{n}^{*}\left( {f_{02},x} \right)}{E\left( {f_{02},x} \right)}{\frac{f_{02}}{f_{e}}}\ {{x\left( {f_{e} < \frac{1}{4\; A}} \right)}}}} + {\int_{- \frac{\pi}{2}}^{\frac{\pi}{2}}{{H_{l}\left( {f_{03},x} \right)}{H_{n}^{*}\left( {f_{03},x} \right)}{E\left( {f_{03},x} \right)}{\frac{f_{03}}{f_{e}}}\ {x}}} + {\int_{- \pi}^{- \frac{\pi}{2}}{{H_{l}\left( {f_{01},x} \right)}{H_{n}^{*}\left( {f_{01},x} \right)}{E\left( {f_{01},x} \right)}{\frac{f_{01}}{f_{e}}}\ {x}}}}} & (10)\end{matrix}$

In Equation (10), the second to fourth terms on the right side representcontribution at a time of a following wave, in other words, representthe degree of the frequency component of the wave when running on thefallowing wave included in the cross spectrum. The parameter A, thethree encounter frequencies f₀₁, f₀₂, and f₀₃ corresponding to theabsolute frequency, and the Yacoubian are each defined as shown inEquation (11) below. Note that, in Equation (11), U is the vessel speedand g is the gravitational acceleration.

$\begin{matrix}\left\lbrack {{Mathematical}\mspace{14mu} {Expression}\mspace{14mu} 11} \right\rbrack & \; \\\left\{ \begin{matrix}{A = {\frac{2\pi}{g}U\; \cos \; x}} \\{f_{01} = \frac{1 - \sqrt{1 - {4\; {Af}_{e}}}}{2\; A}} \\{f_{02} = \frac{1 + \sqrt{1 - {4\; {Af}_{e}}}}{2\; A}} \\{f_{03} = \frac{1 + \sqrt{1 + {4\; {Af}_{e}}}}{2\; A}} \\{{\frac{f_{01}}{f_{e}}} = \frac{1}{\sqrt{1 - {4\; {Af}_{e}}}}} \\{{\frac{f_{02}}{f_{e}}} = \frac{1}{\sqrt{1 - {4\; {Af}_{e}}}}} \\{{\frac{f_{03}}{f_{e}}} = \frac{1}{\sqrt{1 + {4\; {Af}_{e}}}}}\end{matrix} \right. & (11)\end{matrix}$

Here, in the case that the integration range with respect to theencounter angle x is divided into a sufficiently large number K of finesections, the response function and the directional wave spectrum of avariation amount can be constant in each of the small sections.Therefore, Equation (10) can be discretized into Equation (12) below.Note that, in Equation (12), K1 (where, 0≦K1≦K/2) is the number of thefine sections which are in a following wave state in the discreteintegration range.

$\begin{matrix}{\mspace{79mu} \left\lbrack {{Mathematical}\mspace{14mu} {Expression}\mspace{14mu} 12} \right\rbrack} & \; \\{{{\varphi_{\ln}\left( f_{e} \right)} = {{\Delta \; x{\sum\limits_{k = 1}^{K}\; {{H_{l,k}\left( f_{01} \right)}{H_{n,k}^{*}\left( f_{01} \right)}{E_{k}\left( f_{01} \right)}}}} + {\Delta \; x{\sum\limits_{k = 1}^{K\; 1}\; {{H_{l,k}\left( f_{02} \right)}{H_{n,k}^{*}\left( f_{02} \right)}{E_{k}\left( f_{02} \right)}}}} + {\Delta \; x{\sum\limits_{k = 1}^{K\; 1}\; {{H_{l,k}\left( f_{03} \right)}{H_{n,k}^{*}\left( f_{03} \right)}{E_{k}\left( f_{03} \right)}}}}}}\mspace{20mu} {Where}\mspace{20mu} \left\{ \begin{matrix}{{\Delta \; x} = {2{\pi/K}}} \\{{E_{k}\left( f_{0*} \right)} = {E_{k}\left( {f_{0*},x_{k}} \right)}} \\{x_{k} = {{- \pi} + {\left( {k - 1} \right)\Delta \; x}}} \\{{H_{l,k}\left( f_{0*} \right)} = {H_{l}\left( {f_{0*},x_{k}} \right)}} \\{{H_{n,k}^{*}\left( f_{0*} \right)} = {H_{n}^{*}\left( {f_{0*},x_{k}} \right)}}\end{matrix} \right.} & (12)\end{matrix}$

Here, in the case that the pitching angle, the rolling angle, and theheaving displacement are respectively any given variation amounts θ, φ,and η, the cross spectrum Φ(f_(e)) is a 3×3 matrix, and Equation (12)can be expressed in a matrix as Equation (13) below. Note that, inEquation (13), H(f₀₁) is a 3×K matrix, H(f₀₂) and H(f₀₃) are 3×K1matrices, E(f₀₁) is K×K diagonal matrix, E(f₀ 2) and E(f₀₃) are K1×K1diagonal matrices. Further, the symbol T represents a transposed matrix.

$\begin{matrix}{\mspace{79mu} \left\lbrack {{Mathematical}\mspace{14mu} {Expression}\mspace{14mu} 13} \right\rbrack} & \; \\{{{\Phi \left( f_{e} \right)} = {{{H\left( f_{01} \right)}{E\left( f_{01} \right)}{H\left( f_{01} \right)}^{*T}} + {{H\left( f_{02} \right)}{E\left( f_{02} \right)}{H\left( f_{02} \right)}^{*T}} + {{H\left( f_{03} \right)}{E\left( f_{03} \right)}{H\left( f_{03} \right)}^{*T}}}}\mspace{20mu} {Where}\mspace{20mu} {{\Phi \left( f_{e} \right)} = \begin{pmatrix}{\Phi_{\theta\theta}\left( f_{e} \right)} & {\Phi_{\theta\varphi}\left( f_{e} \right)} & {\Phi_{\theta\eta}\left( f_{e} \right)} \\{\Phi_{\varphi\theta}\left( f_{e} \right)} & {\Phi_{\varphi\varphi}\left( f_{e} \right)} & {\Phi_{\varphi\eta}\left( f_{e} \right)} \\{\Phi_{\eta\theta}\left( f_{e} \right)} & {\Phi_{\eta\varphi}\left( f_{e} \right)} & {\Phi_{\eta\eta}\left( f_{e} \right)}\end{pmatrix}}\mspace{20mu} {{{H\left( f_{01} \right)} = \begin{pmatrix}{H_{\theta 1}\left( f_{01} \right)} & \ldots & {H_{\theta \; K}\left( f_{01} \right)} \\{H_{\varphi 1}\left( f_{01} \right)} & \ldots & {H_{\varphi \; K}\left( f_{01} \right)} \\{H_{\eta 1}\left( f_{01} \right)} & \ldots & {H_{\eta \; K}\left( f_{01} \right)}\end{pmatrix}},\mspace{20mu} {{E\left( f_{01} \right)} = \begin{pmatrix}{E_{1}\left( f_{01} \right)} & \ldots & 0 \\\vdots & \ddots & \vdots \\0 & \ldots & {E_{K}\left( f_{01} \right)}\end{pmatrix}}}\mspace{20mu} {{{H\left( f_{0\; i} \right)} = \begin{pmatrix}{H_{\theta 1}\left( f_{0\; i} \right)} & \ldots & {H_{\theta \; K\; 1}\left( f_{0\; i} \right)} \\{H_{\varphi 1}\left( f_{0\; i} \right)} & \ldots & {H_{\varphi \; K\; 1}\left( f_{0\; i} \right)} \\{H_{\eta 1}\left( f_{0\; i} \right)} & \ldots & {H_{\eta \; K\; 1}\left( f_{0\; i} \right)}\end{pmatrix}},\mspace{20mu} {{E\left( f_{0\; i} \right)} = {{\begin{pmatrix}{E_{1}\left( f_{0\; i} \right)} & \ldots & 0 \\\vdots & \ddots & \vdots \\0 & \ldots & {E_{K\; 1}\left( f_{0\; i} \right)}\end{pmatrix}\mspace{31mu} i} = 2}},3}} & (13)\end{matrix}$

Because the cross spectrum matrix Φ(f_(e)) is an Hermitian matrix, it isenough to deal with the upper triangular matrix. Further, in the casethat Equation (13) is expressed by separating the real part and theimaginary part and by introducing an error term W associated withobservation, Equation (13) is expressed by a linear regression modelrepresented by Equation (14) below.

[Mathematical Expression 14]

y=Ax+W   (14)

In Equation (14), y is a vector constituted by the real part and theimaginary part of the cross spectrum matrix Φ(f_(e)). The parameter A isa coefficient matrix constituted by a logical value of the responsefunction of the hull motion in wave. The error term W is a white noisehaving a statistical characteristics of the average 0 and following avariance-covariance matrix Σ. The encounter angle x is an unknown vectorconstituted by the discretized directional wave spectrum.

In Equation (14), it is assumed that the cross spectrum is obtained intime series, the directional wave spectrum can be estimated in timeseries. This corresponds to considering Equation (14) as a time varyingsystem, and Equation (14) can thus be extended into Equation (15) below,where time is represented by index t.

[Mathematical Expression 15]

y _(t) =A _(t) x _(t) +W _(t)   (15)

Equation (15) is formally equivalent to an observation model in ageneral state-space model. Therefore, by introducing as a system model asmoothing prior distribution, in which the directional wave spectrumchanges smoothly with time, (see Equation (16) below), the problem ofestimating the directional wave spectrum can be treated as a problem ofthe state estimation of the general state-space model shown by Equation(16) below.

$\begin{matrix}\left\lbrack {{Mathematical}\mspace{14mu} {Expression}\mspace{14mu} 16} \right\rbrack & \; \\\left\{ \begin{matrix}{x_{t} = {x_{t - 1} + v_{t}}} \\{y_{t} = {{A_{t}x_{t}} + W_{t}}}\end{matrix} \right. & (16)\end{matrix}$

In Equation (16), x_(t) is a state vector, v_(t) is a system noisevector, y_(t) is an observation vector, A_(t) is a state transitionmatrix, and W_(t) is an observation noise vector. Here, considering thatthe directional wave spectrum is not negative, the logarithm of thestate vector x_(t) is replaced anew by x_(t) to deform Equation (16)into a general state-space model represented by Equation (17) below.

$\begin{matrix}\left\lbrack {{Mathematical}\mspace{14mu} {Expression}\mspace{14mu} 17} \right\rbrack & \; \\\left\{ \begin{matrix}{x_{t} = {x_{t - 1} + v_{t}}} \\{y_{t} = {{A_{t}{F\left( x_{t} \right)}} + W_{t}}}\end{matrix} \right. & (17)\end{matrix}$

Here, F(x_(t)) means that F(x_(t)) is exponential to all the elements.Further, the elements of the state vector are configured as Equation(18) below.

$\begin{matrix}\left\lbrack {{Mathematical}\mspace{14mu} {Expression}\mspace{14mu} 18} \right\rbrack & \; \\\left\{ \begin{matrix}{{F\left( x_{t} \right)} = {\exp \left\lbrack x_{t} \right\rbrack}} \\{x_{l}^{T} = \left\lbrack {{\ln \left( {x_{1},t} \right)},\ldots \mspace{14mu},{\ln \left( {x_{J},t} \right)}} \right\rbrack} \\{{F\left( x_{t} \right)}^{T} = {\exp \left\lbrack {{\ln \left( {x_{1},t} \right)},\ldots \mspace{14mu},{\ln \left( {x_{J},t} \right)}} \right\rbrack}} \\{{x_{j,t} = {E_{k,t}\left( f_{0\; i} \right)}},\mspace{14mu} {i = {1 - m}},\mspace{14mu} {j = {1 - J}},\mspace{14mu} {J = {m \times k}},\mspace{14mu} {k = {1 - {K.}}}}\end{matrix} \right. & (18)\end{matrix}$

In Equation (18), m is the number of division of the absolute frequency.Equation (17) is a non-linear observation model, in other words, anon-linear state space model. Therefore, in order to estimate the state,it is necessary to use a method effective in non-linear filtering.Conventionally, a particle filter is used, but this filter has a verylarge calculation load. To address this issue, a state estimation methodby using an Ensemble Kalman Filter is introduced in the presentembodiment; however, an Ensemble Kalman Filter cannot be applied toEquation (17) in a non-linear observation model as it is. To solve thisproblem, the extended state vector represented by Equation (19) below isconsidered.

$\begin{matrix}\left\lbrack {{Mathematical}\mspace{14mu} {Expression}\mspace{14mu} 19} \right\rbrack & \; \\{Z_{i} = \begin{pmatrix}x_{t} \\{A_{t}{F\left( x_{t} \right)}}\end{pmatrix}} & (19)\end{matrix}$

Further, the extended observation matrix and the extended statetransition matrix represented by Equation (20) below are considered.

$\begin{matrix}\left\lbrack {{Mathematical}\mspace{14mu} {Expression}\mspace{14mu} 20} \right\rbrack & \; \\\left\{ \begin{matrix}{{{\overset{\sim}{A}}_{t} = \left( {O_{l \times k}\mspace{14mu} I_{l \times l}} \right)},} \\{{{\overset{\sim}{f}}_{t}\left( {z_{t - 1},v_{t}} \right)} = \begin{pmatrix}{x_{t - 1} + v_{t}} \\{A_{t}{F\left( {x_{i - 1} + v_{i}} \right)}}\end{pmatrix}}\end{matrix} \right. & (20)\end{matrix}$

As a result, Equation (21) below holds for x_(t), and an extended systemmodel is obtained. Further, regarding to y_(t), a formally linearextended observation model represented by Equation (22) below can beobtained. Because x_(t) and y_(t) are an extended state space model inthe linear observation, the state estimation by an Ensemble KalmanFilter can be realized. Note that the application of the Ensemble KalmanFilter is a known technique and is thus not described here.

[Mathematical Expression 21]

z _(t) =f _(t)(z _(t-1) , v _(i))   (21)

[Mathematical Expression 22]

y _(t) =Ā _(t) z _(t) +w _(t)   (22)

By the above-described process of step S7, the sea phenomena estimationsection 24 probability statistically calculates the directional wavespectrum on the basis of the hull response calculated in step S5 and thecross spectra of the respective hull motions calculated in step S6 (stepS7).

In the process of step S7, by probability statistically processing thehull response in a predetermined period from past to now and thetime-series data of the cross spectra of the respective hull motions,the directional wave spectrum is calculated in real time. Thus, thedirectional wave spectrum with high accuracy can be derived.

Further, by this method, the directional wave spectrum is estimated onthe basis of the state estimation by using the Ensemble Kalman Filter;thus, it is possible to estimate the directional wave spectrum with highaccuracy in much shorter calculation time than in the method using theconventional Monte Carlo filter.

Note that when the process of step S7 has finished, the process goes tostep S8, and the sea phenomena estimation section 24 estimates the seaconditions on the basis of the directional wave spectrum calculated instep S7 (step S8). In step S8, it is possible to estimate the seaconditions such as the wave direction, the wave cycle, the significantwave height, and the like in the local water in which the vessel isrunning, on the basis of the directional wave spectrum calculated instep S7.

An embodiment of the present invention is described above, but the aboveembodiment describes merely one of the application examples of thepresent invention. Therefore, it is not intended that the technicalscope of the present invention is limited to the specific configurationof the above embodiment.

REFERENCE SIGNS LIST

-   -   1: Sea phenomena estimation system    -   2: Satellite compass    -   3: Information processing device (Transverse metacenter height        estimation device)    -   4: Display    -   21: Measurement section    -   22: Time history memory    -   23: Hull state data calculation section (Transverse metacenter        height estimation section)    -   24: Sea phenomena estimation section    -   25: Output section

1. A transverse metacenter height estimation device, comprising: a timehistory memory unit that stores time-series data of a rolling angle of ahull; and a transverse metacenter height estimation unit that estimatesa transverse metacenter height of the hull based on the time-series dataof the rolling angle of the hull stored in the time history memory unit,wherein the transverse metacenter height estimation unit calculates arolling natural frequency based on the time-series data of the rollingangle of the hull, and estimates, while using the calculated rollingnatural frequency as an observation model, the transverse metacenterheight by performing state estimation based on a general state-spacemodel in which the transverse metacenter height and a radius of gyrationof the hull are state variables.
 2. The transverse metacenter heightestimation device according to claim 1, wherein the transversemetacenter height estimation unit estimates the transverse metacenterheight GM by assuming a transverse metacenter height GM_(n) at time nand a radius of gyration kn at time n as state variables in Equation (1)representing a relationship between the rolling natural frequency andthe transverse metacenter height: $\begin{matrix}\left\lbrack {{Mathematical}\mspace{14mu} {Expression}\mspace{14mu} 1} \right\rbrack & \; \\{T = {\frac{1}{f} = \frac{2\pi \; k}{\sqrt{gGM}}}} & (1)\end{matrix}$ where T is a rolling natural period, f is a rollingnatural frequency, k is a radius of gyration, g is a gravitationalacceleration, and GM is a transverse metacenter height, and byperforming state estimation based on a general state-space model whichis represented by Equation (2): $\begin{matrix}\left\lbrack {{Mathematical}\mspace{14mu} {Expression}\mspace{14mu} 2} \right\rbrack & \; \\\left\{ {\begin{matrix}{x_{n} = {x_{n - 1} + v_{n}}} \\{T_{n} = {h\left( {x_{n},w_{n}} \right)}}\end{matrix}{Where}\left\{ \begin{matrix}{x_{n} = \left( {{GM}_{n}\mspace{14mu} k_{n}} \right)^{T}} \\{v_{n} = \left( {v_{1,n}\mspace{14mu} v_{2,n}} \right)^{T}}\end{matrix} \right.} \right. & (2)\end{matrix}$ where Vn is a system noise, Wn is an observation noise,and T_(n) is a rolling natural period and is observed data.
 3. Atransverse metacenter height estimation method comprising the steps of;storing time-series data of a rolling angle of a hull in a time historymemory; and estimating a transverse metacenter height of the hull basedon the time-series data of the rolling angle of the hull stored in thetime history memory, wherein in the step of estimating the transversemetacenter height of the hull, a rolling natural frequency is calculatedbased on the time-series data of the rolling angle of the hull, and thetransverse metacenter height is estimated, while the calculated rollingnatural frequency is used as an observation model, by performing stateestimation based on a general state-space model in which the transversemetacenter height and an radius of gyration of the hull are statevariables.